Starry Night?

An exploration of nighttime lights and the Dark Sky Movement in the US

Author

Carly Staebell

Published

December 22, 2025

Introduction

Electric lights have benefited society immensely since their invention in 1879. Modern work, recreation, and transportation are made possible largely due to artificial lighting. However, when artificial outdoor lighting becomes inefficient and unnecessary, it is known as light pollution, with negative impacts on both wildlife and humans [1].

Concerns about the impacts of light pollution on astronomical observation led to the 1988 formation of the nonprofit now known as DarkSky International [2]. Since then, DarkSky has grown into a recognized global authority on light pollution with a mission to safeguard communities and wildlife against light pollution [3]. The organization’s International Dark Sky Places program allows communities, parks, and protected areas to seek certification for efforts to preserve and protect their dark skies [4]. While many communities implement outdoor lighting policies without seeking recognition, a DarkSky certification is an indicator of places that have made dark sky conservation a priority.

This project aims to explore trends in nighttime light levels first from a national contiguous US perspective, and through more detailed case studies of individual cities. Three DarkSky certified communities [4] with different sizes and designation dates have been selected:

Flagstaff, Arizona: Population ~77,000, designated 2001;
Homer Glen, Illinois: Population ~25,000, designated 2011;
Bee Cave, Texas: Population ~9,000, designated 2023.

These three communities were each paired with a community of similar size that has fewer, if any, ordinances controlling light pollution. Each community pair was examined to see if there are observable differences in nighttime radiance.

Materials and methods

The following data were used.

NASA’s Black Marble annual nighttime lights (VNP46A4 data product) [5] sourced via the blackmarbler package [6]
- Spatial resolution: 500 m
US Census Bureau TIGER data [7] sourced via the tigris package [8]
US Census Bureau American Community Survey (ACS) population data [9] sourced via the tidycensus package [10]
MODIS Type 1 yearly land use/land cover (MCD12Q1 data product) sourced via AppEEARS [11]
- Spatial resolution: 500 m

Code

# Load libraries
library(tidyverse)
library(leaflet)
library(kableExtra)
library(htmlwidgets)
library(widgetframe)
knitr::opts_chunk$set(widgetframe_widgets_dir = 'widgets' ) 
knitr::opts_chunk$set(cache=TRUE)  # cache the results for quick compiling

library(blackmarbler)
library(sf)
library(terra)
library(tigris)
library(tidyterra)
library(tidycensus)
library(patchwork)
library(piggyback)
library(ggiraph)
library(scales)

Download and clean all required data

National and statewide nighttime lights

A US Census TIGER shapefile was downloaded to provide geographic state boundaries. A separate script was used to download and save annual nighttime lights rasters for the contiguous US (‘black_marble_data_download.R’). The script was also used to obtain extracted nighttime lights values aggregated at the national and state levels. Nighttime lights are provided as radiance values with units of $\frac{n W}{cm^2 sr}$.

Code

# Download and process CONUS NTL data
# Define national spatial extent
nonconus <- c("Guam", "Hawaii", "Alaska",
              "Commonwealth of the Northern Mariana Islands",
              "United States Virgin Islands", "American Samoa", "Puerto Rico")

states_sf <- states() |>
  filter(!NAME %in% nonconus) |>
  select(NAME, geometry) |>
  st_transform(crs = "epsg:4326")

# Dissolve individual state boundaries to get CONUS outline
conus <- st_union(states_sf)|>
  st_as_sf()

## Read in national raster data for 2024 (to be used for leaflet map later)
pb_download("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif", 
            repo = "cstaebell/final-project-cstaebell")

us_2024_rast <- rast("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif") |>
  mask(conus)

# create log transformed version for plotting
log_us_2024 <- app(us_2024_rast, fun = log1p)

# Downsample log data for map (commented out since file written)
#aggregate(log_us_2024, fact = 4, fun = "median", filename = "data/log_us_2024_downsampled.tif", overwrite = TRUE)

#---
# Download and process annual US mean NTL values
#---
years = 2014:2024
conus_means = numeric(length(years))

for (i in 1:length(years)) {
  pb_download(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
                             years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell")
  annual_mean <- readRDS(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
                             years[i], ".Rds", sep = "")) |>
    pull(ntl_mean)
  
  conus_means[i] <- annual_mean
}

conus_df <- data.frame(year = years, mean_ntl = conus_means)

#---
# Download and process mean NTL values for states
#---
state_means <- matrix(nrow = 49, ncol = length(years))

# Read in annual data for states and create a matrix of state means
for (i in 1:length(years)) {
  pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
                             years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell")
  annual_mean <- readRDS(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
                             years[i], ".Rds", sep = ""))
  
  state_means[,i] <- annual_mean[,2]

  # check that state names are ordered the same way for each year
  if (i == 1) {
    state_names <- annual_mean[,1]
  } else {
    if (sum(state_names != annual_mean[,1]) != 0) {
      print("Different state names")
    }
  }
}

# Create data frame of state means and add percent change column
state_means_df <- data.frame(state_means)
colnames(state_means_df) <- paste("ntl_", years, sep = "")
state_means_df <- state_means_df |>
  mutate(state = state_names, 
         pct_change = (ntl_2024 - ntl_2014) / ntl_2014 * 100)


# Combine state NTL data with geographic data
states_ntl <- states_sf |>
  inner_join(state_means_df, by = c("NAME" = "state"))

Case Studies

For case study analyses, population data and local geographies for each Dark Sky Place were downloaded from the US Census.

Code

# Download population data
# US population to search for comparison cities
us_pop <- get_acs(geography = "place", state = NULL, 
                  variables = "B01003_001", key = census_api_key)

# Download local geographies: selected Dark Sky Places
flagstaff_sf <- places(state = "AZ", cb = TRUE, year = 2024) |>
  filter(NAME == "Flagstaff") |>
  st_transform(crs = "epsg:4326") |>
  left_join(us_pop, by = "GEOID")

homerglen_sf <- places(state = "IL", cb = TRUE, year = 2024) |>
  filter(NAME == "Homer Glen") |>
  st_transform(crs = "epsg:4326") |>
  left_join(us_pop, by = "GEOID")

beecave_sf <- places(state = "TX", cb = TRUE, year = 2024) |>
  filter(NAME == "Bee Cave") |>
  st_transform(crs = "epsg:4326") |>
  left_join(us_pop, by = "GEOID")

ACS population data were filtered to return a listing of candidate municipalities with populations within 1% of each Dark Sky community population. Final Comparison cities were selected based on population match and state, with preference given to cities within states that are not known to have laws regarding outdoor lighting geared toward dark sky protection [12].

Code

# Get listing of suitable comparison cities
flagstaff_sim_pop <- us_pop |> filter(estimate >= (flagstaff_sf$estimate - 0.005*flagstaff_sf$estimate) 
                                & estimate <= (flagstaff_sf$estimate + 0.005*flagstaff_sf$estimate)) |>
  separate(NAME, into = c("city", "state"), sep = ", ")

homerglen_sim_pop <- us_pop |> filter(estimate >= (homerglen_sf$estimate - 0.005*homerglen_sf$estimate) 
                                & estimate <= (homerglen_sf$estimate + 0.005*homerglen_sf$estimate)) |>
  separate(NAME, into = c("city", "state"), sep = ", ")

beecave_sim_pop <- us_pop |> filter(estimate >= (beecave_sf$estimate - 0.005*beecave_sf$estimate) 
                                & estimate <= (beecave_sf$estimate + 0.005*beecave_sf$estimate)) |>
  separate(NAME, into = c("city", "state"), sep = ", ")

# Download/process comparison city data
 # Flagstaff counterpart
scranton_sf <- places(state = "PA", cb = TRUE, year = 2024) |>
  filter(NAME == "Scranton") |>
  left_join(us_pop, by = "GEOID") |>
  st_transform(crs = "epsg:4326")

 # Homer Glen counterpart
belvidere_sf <- places(state = "IL", cb = TRUE, year = 2024) |>
  filter(NAME == "Belvidere") |>
  left_join(us_pop, by = "GEOID") |>
  st_transform(crs = "epsg:4326")

 # Bee Cave counterpart
clanton_sf <- places(state = "AL", cb = TRUE, year = 2024) |>
  filter(NAME == "Clanton") |>
  left_join(us_pop, by = "GEOID") |>
  st_transform(crs = "epsg:4326")

city_names <- c("flagstaff", "homerglen", "beecave", "scranton", 
                           "belvidere", "clanton")

The final city pairs are as follows:

Pair	Community	Population (2023)
A	Flagstaff, AZ	76,333
A	Scranton, PA	76,074
B	Homer Glen, IL	24,516
B	Belvidere, IL	24,510
C	Bee Cave, TX	8,861
C	Clanton, AL	8,862

Annual population data for each city were downloaded from the ACS. The 5-year ACS was used for years 2014-2023. Since the 1-year ACS is only available for cities with populations > 65,000, 2024 data could only be obtained for Pair A cities. For the remaining cities, 2023 data was used for 2024.

Code

# 2014-2023 uses 5-year ACS population data
annual_acs <- purrr::map(years[1:10],
               ~ get_acs(geography = "place", 
                         state = c("AZ", "IL", "TX", "PA", "AL"),
                         variables = "B01003_001", 
                         key = census_api_key, 
                         year = .x, 
                         survey = "acs5")) |>
  map2(years[1:10], ~ mutate(.x, id = .y)) |>
  list_rbind() |>
  filter(GEOID %in% c(flagstaff_sf$GEOID, homerglen_sf$GEOID, beecave_sf$GEOID,
                      scranton_sf$GEOID, belvidere_sf$GEOID, clanton_sf$GEOID)) |>
  separate(NAME, into = c("city", "state"), sep = ", ") |>
  rename(year = id, pop_est = estimate) |>
  mutate(city = case_when(
    city == "Flagstaff city" ~ "flagstaff",
    city == "Scranton city" ~ "scranton",
    city == "Homer Glen village" ~ "homerglen",
    city == "Belvidere city" ~ "belvidere",
    city == "Bee Cave city" ~ "beecave",
    city == "Clanton city" ~ "clanton",
  ))

# 2024 requires 1-year ACS: only available for geographies with populations of 65000 and greater
pop_2024 <- get_acs(geography = "place", state = c("AZ", "PA"),
                                variables = "B01003_001", key = census_api_key, 
                                year = 2024, survey = "acs1")

flagstaff_2024 <- pop_2024 |>
  filter(GEOID == flagstaff_sf$GEOID) |>
  separate(NAME, into = c("city", "state"), sep = ", ") |>
  rename(pop_est = estimate) |>
  mutate(year = 2024, city = "flagstaff")

scranton_2024 <- pop_2024 |>
  filter(GEOID == scranton_sf$GEOID) |>
  separate(NAME, into = c("city", "state"), sep = ", ") |>
  rename(pop_est = estimate) |>
  mutate(year = 2024, city = "scranton")

# Since data unavailable for smaller cities for 2024, impute 2023 data
other_cities_2024 <- annual_acs |>
  filter(year == 2023 & city != "flagstaff" & city != "scranton") |>
  mutate(year = 2024)

all_annual_acs <- annual_acs |>
  rbind(flagstaff_2024, scranton_2024, other_cities_2024)

After selecting cities, the annual US nighttime lights rasters were read in and subsetted to create a raster and dataframe of mean values for each city. Since simply looking at annual means for each city may be more of a proxy for overall development than for lighting policy, land use/land cover data¹ were also brought in to help make comparisons across communities more meaningful. The land cover rasters were cropped to each city, and then the urban areas were isolated, allowing for calculation of annual mean nighttime lights from urban areas only.

Code

# Nighttime lights and land cover data are in a single chunk to avoid external pointer errors when rendering

#---
# Nighttime lights
#---
# Subset nighttime lights rasters to each city
city_sf_list <- list(flagstaff_sf, homerglen_sf, beecave_sf, 
                     scranton_sf, belvidere_sf, clanton_sf)
city_ntl_rasters <- vector("list", length = length(city_names))

# Load files for each year
for (i in 1:length(years)){
  pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t", 
                     years[i], ".tif", sep = ""), repo = "cstaebell/final-project-cstaebell")
  rast_file <- paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t", 
                     years[i], ".tif", sep = "")
  
  year_rast <- rast(rast_file)
  
  # Loop over cities
  for (j in 1:length(city_sf_list)) {
    # Crop and mask to city extent
    current_rast <- year_rast |>
      terra::crop(city_sf_list[[j]]) |>
      mask(city_sf_list[[j]])
    
    # Create list entry during first iteration, then append in subsequent iterations
    if (i == 1) {
    city_ntl_rasters[[j]] <- list(current_rast)
    } else {
    city_ntl_rasters[[j]][[length(city_ntl_rasters[[j]]) + 1]] <- current_rast
    }
  }
}

# Stack city lists 
city_ntl_rasters <- lapply(city_ntl_rasters, rast)

# Create separate rasters for each city 
#city_ntl_var_names <- paste0(city_names, "_ntl_rast")
#for (c in 1:length(city_ntl_var_names)) {
#  assign(city_ntl_var_names[c], city_ntl_rasters[[c]])
#}

# Calculate annual means
annual_ntl_means <- vector("list", length = length(city_names))
for (i in 1:length(city_ntl_rasters)) {
  ntl_vals <- data.frame(values(city_ntl_rasters[[i]])) |>
    summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
    pivot_longer(cols = t2014:t2024) |>
    mutate(name = years, city = city_names[i]) |>
    rename(year = name, mean_ntl = value)
  
  annual_ntl_means[[i]] <- ntl_vals
}

ntl_means_df <- list_rbind(annual_ntl_means) |>
  left_join(all_annual_acs, by = join_by(city == city, year == year)) |>
  mutate(percap_mean_ntl = mean_ntl / pop_est)

#---
# Land cover data
#---
# Download land cover data 
pb_download("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif", 
            repo = "cstaebell/final-project-cstaebell")
pb_download("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif", 
            repo = "cstaebell/final-project-cstaebell")
land_use_rast_2019 <- rast("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif")
land_use_rast_2020 <- rast("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif")

# Subset land cover rasters to each city
land_use_rasters <- vector("list", length = length(city_sf_list))
urban_ntl_dfs <- vector("list", length = length(city_sf_list))
urban_ntl_vals <- vector("list", length = length(city_sf_list))
group_assign <- c("Dark Sky", "Dark Sky", "Dark Sky", "Comparison", "Comparison", "Comparison")
pair_assign <- c("A", "B", "C", "A", "B", "C")


for (i in 1:length(city_names)) {
  # Specify which data to use (2019 for all cities but Clanton)
  if (i == 6) {
    land_use_rast <- land_use_rast_2020
  } else {
    land_use_rast <- land_use_rast_2019
  }
  
  # Crop and mask to city extent
  current_rast <- land_use_rast |>
    terra::crop(city_sf_list[[i]]) |>
    mask(city_sf_list[[i]])
  # Save to list for future use 
  land_use_rasters[[i]] <- current_rast
  
  # Mask all but urban land use (land use code = 13)
  urban_area <- current_rast
  urban_area[urban_area != 13] <- NA
  city_urban <- city_ntl_rasters[[i]] |>
    mask(urban_area) 
  
  # Calculate mean NTL for urban areas 
  urban_ntl_dfs[[i]] <- data.frame(values(city_urban)) |>
    summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
    pivot_longer(cols = t2014:t2024) |>
    mutate(name = years, city = city_names[i]) |>
    rename(year = name, mean_urban_ntl = value) 
  
  # Compile values for statistical tests
  urban_ntl_vals[[i]] <- data.frame(values(city_urban)) |>
    pivot_longer(cols = t2014:t2024) |>
    filter(!is.na(value)) |>
    mutate(city = city_names[[i]], group = group_assign[i], pair = pair_assign[i])
}

# Combine all mean NTL data into single data frame
all_ntl_means_df <- ntl_means_df |>
  left_join(list_rbind(urban_ntl_dfs), join_by(year == year, city == city)) |>
  mutate(mean_percap_urban_ntl = mean_urban_ntl / pop_est, 
         pair = case_when(
           city == "flagstaff" | city == "scranton" ~ "A",
           city == "homerglen" | city == "belvidere" ~ "B",
           city == "beecave" | city == "clanton" ~ "C"),
         group = case_when(
           city == "flagstaff" | city == "homerglen" | city == "beecave" ~ "Dark Sky",
           city == "scranton" | city == "belvidere" | city == "clanton" ~ "Comparison")
         )

# Prepare data for sample land cover plots
# Define land cover encoding
Land_Cover_Type_1 <- c(
  "Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2,
  "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5,
  "Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9,
  "Grasslands" = 10, "Permanent wetlands" = 11, "Croplands" = 12, "Urban & built-up" = 13,
  "Cropland/Natural vegetation mosaic" = 14, "Snow & ice" = 15, "Barren/Sparsely vegetated" = 16
)

lcd <- data.frame(
  ID = Land_Cover_Type_1,
  landcover = names(Land_Cover_Type_1),
  stringsAsFactors = FALSE)

# Create data frames for plotting
flagstaff_lulc <- as.data.frame(land_use_rasters[[1]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
  mutate(city = "Flagstaff")

scranton_lulc <- as.data.frame(land_use_rasters[[4]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
  mutate(city = "Scranton")

homerglen_lulc <- as.data.frame(land_use_rasters[[2]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
  mutate(city = "Homer Glen")

belvidere_lulc <- as.data.frame(land_use_rasters[[5]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
  mutate(city = "Belvidere")

beecave_lulc <- as.data.frame(land_use_rasters[[3]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
  mutate(city = "Bee Cave")

clanton_lulc <- as.data.frame(land_use_rasters[[6]], xy = TRUE) |>
  left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2020001000000_aid0001 == ID)) |>
  rename(lc_code = MCD12Q1.061_LC_Type1_doy2020001000000_aid0001) |>
  mutate(city = "Clanton")

all_lulc <- bind_rows(flagstaff_lulc, scranton_lulc, homerglen_lulc, belvidere_lulc, 
                      beecave_lulc, clanton_lulc)
all_lulc$city <- factor(all_lulc$city, levels = c("Flagstaff", "Scranton",
                                                     "Homer Glen", "Belvidere",
                                                     "Bee Cave", "Clanton"))

Statistical testing was done on the nonaggregated urban lights values. Although there are thousands of samples, radiance data exhibit a heavy right skew. The Shapiro-Wilk test was used to confirm non-normality, and since the data are not normal, a Wilcoxon test was used for each community pair to test if the Comparison group radiance is greater than that of the Dark Sky group.

Code

# Prep values for statistical testing
pair_a_uvals <- bind_rows(urban_ntl_vals[[1]], urban_ntl_vals[[4]])
pair_b_uvals <- bind_rows(urban_ntl_vals[[2]], urban_ntl_vals[[5]])
pair_c_uvals <- bind_rows(urban_ntl_vals[[3]], urban_ntl_vals[[6]])

for (i in 1:length(urban_ntl_vals)) {
  shapiro_test <- shapiro.test(urban_ntl_vals[[i]]$value)
  if (shapiro_test$p.value < 0.05) {
    print("Reject normality")
  } else {
    print("Fail to reject normality")
  }
}

wilcox_a <- wilcox.test(value ~ group, data = pair_a_uvals, alternative = "greater")

wilcox_b <- wilcox.test(value ~ group, data = pair_b_uvals, alternative = "greater")

wilcox_c <- wilcox.test(value ~ group, data = pair_c_uvals, alternative = "greater")

# Repeat analysis for only 2021 - 2024 due to population trends
pair_a_recent <- pair_a_uvals |>
  filter(name %in% c("t2021", "t2022", "t2023", "t2024"))
pair_b_recent <- pair_b_uvals |>
  filter(name %in% c("t2021", "t2022", "t2023", "t2024"))
pair_c_recent <- pair_c_uvals |>
  filter(name %in% c("t2021", "t2022", "t2023", "t2024"))

wilcox_a_recent <- wilcox.test(value ~ group, data = pair_a_recent, alternative = "greater")

wilcox_b_recent <- wilcox.test(value ~ group, data = pair_b_recent, alternative = "greater")

wilcox_c_recent <- wilcox.test(value ~ group, data = pair_c_recent, alternative = "greater")

Results

National Trends

Mean nighttime light radiance for the entire contiguous US is plotted below, where a generally increasing trend can be observed over the past ten years.

Code

ggplot(conus_df) +
  geom_point(aes(x = year, y = mean_ntl)) +
  geom_smooth(aes(x = year, y = mean_ntl)) +
  labs(x = "Year", y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
       title = "Annual US Mean Radiance (2014 - 2024)") +
  theme_light()

State means for 2024 are shown in the figure below. Annual radiance tends to be higher for eastern states, which reflects known patterns of population density.

Code

# Create ggplot map for ggiraph
st_2024_mean_map <- ggplot() +
  geom_sf_interactive(data = filter(states_ntl, NAME != "District of Columbia"), 
                      aes(fill = ntl_2024,
                      tooltip = paste(NAME, "\n", round(ntl_2024, 2)),
                      data_id = NAME)) +
  theme_void() +
  theme(text = element_text(family = "Arial")) +
  scale_fill_viridis_c() +
  labs(title = "Mean Nighttime Lights by State, 2024", 
       fill = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
       caption = "Data: NASA Black Marble & US Census Bureau")

# Transform into interactive ggiraph object
girafe(ggobj = st_2024_mean_map,
       options = list(opts_hover(css = "stroke:cyan;stroke-width:2px"),
                      opts_hover_inv(css = "opacity:0.8;"),
                      opts_tooltip(css = "background-color:white;color:black;padding:15px;
                                          border-radius:10px;
                                          font-family:Arial;font-size:15px;"
                                          )))

In the period from 2014 to 2024, most states experienced an increase in detected nighttime radiance. Texas, Montana, Idaho, and Utah all show particularly large increases upwards of 40%. More research is required to explain these trends, but possible drivers include recent population growth and urban development. On the other hand, several states have stayed fairly constant or even decreased over the 10-year span.

Code

st_pct_change_map <- ggplot() + 
  geom_sf_interactive(data = filter(states_ntl, NAME != "District of Columbia"), 
                      aes(fill = pct_change,
                          tooltip = paste(NAME, "\n", round(pct_change, 1),"%", 
                                          "\n2014: ", round(ntl_2014, 2),
                                          "\n2024: ", round(ntl_2024, 2)),
                          data_id = NAME)) +
  theme_void() +
  theme(text = element_text(family = "Arial")) +
  scale_fill_gradient2(low = "navy", 
                       mid = "white",
                       high = "darkred",
                      midpoint = 0) +
  labs(title = "Mean Radiance Percent Change (2014 - 2024)",
       fill = "% Change",
       caption = "Data: NASA Black Marble & US Census Bureau")

girafe(ggobj = st_pct_change_map,
       options = list(opts_hover(css = "stroke:black;stroke-width:2px"),
                 opts_hover_inv(css = "opacity:0.8;"),
                 opts_tooltip(css = "background-color:white;color:black;padding:15px;
                                     border-radius:10px;
                                     font-family:Arial;font-size:15px;")))

Case Studies

Code

# Leaflet prep
fs_coord <- st_centroid(flagstaff_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Flagstaff", pair = "A", group = "Dark Sky")

sc_coord <- st_centroid(scranton_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Scranton", pair = "A", group = "Comparison")

hg_coord <- st_centroid(homerglen_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Homer Glen", pair = "B", group = "Dark Sky")

bv_coord <- st_centroid(belvidere_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Belvidere", pair = "B", group = "Comparison")

bc_coord <- st_centroid(beecave_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Bee Cave", pair = "C", group = "Dark Sky")

cl_coord <- st_centroid(clanton_sf) |>
  st_coordinates() |>
  data.frame() |>
  mutate(city = "Clanton", pair = "C", group = "Comparison")

markers_df <- bind_rows(fs_coord, sc_coord, hg_coord, bv_coord, bc_coord, cl_coord)
markers_df$pair <- as.factor(markers_df$pair)

log_us_2024_ds <- rast("data/log_us_2024_downsampled.tif")

vals <- values(log_us_2024_ds, mat = FALSE)
midpoint <- 3.5

breaks <- c(min(vals, na.rm = TRUE),
            seq(min(vals, na.rm = TRUE), midpoint, length.out = 5)[-1],
            seq(midpoint, max(vals, na.rm = TRUE), length.out = 5)[-1])

pal <- colorBin(palette = c("black", "yellow", "red"),
                domain = vals,
                bins = breaks,
                na.color = "transparent")


city_color_pal <- colorFactor(palette = c("firebrick", "royalblue3", "green4"),
                           domain = markers_df$pair)

city_popup <- paste0("<b>", markers_df$city, "</b><br/>",
                     "<i>", markers_df$group, " Community</i><br/>",
                     "Pair: ", markers_df$pair)

# Leaflet map
leaflet() |>
  addProviderTiles("CartoDB") |>
  addRasterImage(log_us_2024_ds, group = "US NTL", colors = pal) |>
  addCircleMarkers(data = markers_df, ~X, ~Y, fillColor = ~city_color_pal(pair), 
                   fillOpacity = 0.9, color = "black", 
                   label = ~as.character(city), popup = city_popup)

From the plot below, four of the six case study cities exhibit an increasing mean radiance trend similar to that of the national mean. However, the Pair A cities, Flagstaff and Scranton, show decreasing tendencies. The large difference in radiance values within pairs is another notable feature of this plot.

Code

# Prep data frame so that plots show in desired order
all_ntl_means_df$pair <- as.factor(all_ntl_means_df$pair)
pair_names <- c("A" = "Flagstaff & Scranton", "B" = "Homer Glen & Belvidere", "C" = "Bee Cave & Clanton")

# Plot of annual means
ggplot(data = all_ntl_means_df) +
  geom_line(aes(x = year, y = mean_ntl, color = pair, linetype = group), linewidth = 1) +
  facet_wrap(~ pair, labeller = labeller(pair = pair_names)) +
  scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) +
  scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) +
  labs(x = "Year",
       y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
       title = "Mean Annual Radiance",
       linetype = "City Type",
       color = "Pair") +
  theme_light() +
  theme(strip.text.x = element_text(color = "black"),
        plot.title = element_text(hjust = 0.5),
        aspect.ratio = 1.8)

Simply looking at the annual means provides interesting information, but fails to take into account two factors that may influence a city’s observed nighttime light emissions: population over time and land cover. Annual population data reveals that Flagstaff and Bee Cave have both experienced steep population growth from 2014 to 2024. While Belvidere shows a generally decreasing trend, both Belvidere and Homer Glen had relatively stable populations with less than 5% change. Scranton and Clanton’s populations have also remained quite stable.

Code

ggplot(data = all_ntl_means_df) +
  geom_line(aes(x = year, y = pop_est, color = pair, linetype = group), linewidth = 1) +
  facet_wrap(~ pair, labeller = labeller(pair = pair_names), scale = "free_y") +
  scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) +
  scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) +
  labs(x = "Year",
       y = "Population",
       title = "Annual Populations",
       linetype = "City Type",
       color = "Pair") +
  theme_light() +
  theme(strip.text.x = element_text(color = "black"),
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_text(size = 8),
        aspect.ratio = 1.8)

Land cover differences between paired cities helps to explain some of the features and differences observed in the plot of mean nighttime lights versus time. Flagstaff contains much more undeveloped land like grassland and savannas, whereas most of Scranton is classified as urban and built up, which likely accounts for the much lower mean values for Flagstaff. Similarly, Clanton has a very small amount of land classified as urban & built-up, which probably contributes to a low mean radiance.

Code

ggplot(all_lulc, aes(x = x, y = y)) +
  geom_tile(aes(fill = landcover)) +
  facet_wrap(~ city, nrow = 3, scales = "free") +
  coord_quickmap() +
  theme_light() +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        strip.text.x = element_text(color = "black")) +
  labs(fill = "Land Cover",
       caption = "Data: NASA & US Census Bureau") +
  scale_fill_viridis_d()

Given the imbalances in land cover between different cities, it is useful for the sake of comparison to restrict analysis to urban areas. Proceeding this way introduces the assumption that any lighting policy effects outside of non-built-up areas are dwarfed by those on urban areas. From the plots of urban mean radiance, we can see that the radiance for the Dark Sky Places is generally lower than for their paired cities.

Code

# Urban mean NTL plot
ggplot(data = all_ntl_means_df) +
  geom_line(aes(x = year, y = mean_urban_ntl, color = pair, linetype = group), linewidth = 1) +
  facet_wrap(~ pair, labeller = labeller(pair = pair_names), scales = "free_y") +
  scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) +
  scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) +
  labs(x = "Year",
       y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
       title = "Mean Annual Radiance from Urban Areas",
       linetype = "City Type",
       color = "Pair") +
  theme_light() +
  theme(strip.text.x = element_text(color = "black"),
        plot.title = element_text(hjust = 0.5),
        aspect.ratio = 1.8)

The results of the Wilcoxon test on the urban radiance values suggest that the observed differences between the Dark Sky communities and their Comparison counterparts are statistically significant. Given the population discrepancies between the Pair A and Pair C cities in the earlier analysis years, analysis was done on the entire set of data and on a subset of data from 2021-2024, when the populations became more comparable.

Pair	p-value (2014-2024)	p-value (2021-2024)
A	<0.0001	<0.0001
B	<0.0001	<0.0001
C	0.00591	0.03937

Discussion of lighting policies

Pair A

As the home of two major observatories, Flagstaff has been at the forefront of outdoor lighting policy since 1958. The city’s lighting ordinance aims to encourage lighting practices and systems that will minimize light pollution and trespass while still maintaining nighttime safety, security, and productivity. The code includes comprehensive standards for maximum total outdoor light output, shielding requirements, time limits, sign illumination, and lamp types. [13]

Scranton’s Zoning Ordinance includes one short section devoted to lighting, which mainly focuses on curbing light trespass and nuisance glare. Maximum illumination levels are defined for two lighting levels. [14]

Pair B

Homer Glen’s outdoor lighting ordinance is rooted in the desire to preserve rural character and aesthetic value. The ordinance establishes maximum light output and inclination for various outdoor lighting applications, permitted hours for lighting, shielding requirements, and bulb types. [15]

Similar to Scranton, Belvidere includes exterior lighting standards in their code of ordinances but does so for the purposes of promoting traffic safety and preventing nuisances. The ordinance sets limits on illumination intensity and hours of illumination outside of buildings, plus some requirements for light orientation and shielding. However, these provisions are not as strict or clearly defined as those in the Homer Glen ordinanace. [16]

Pair C

Bee Cave’s Unified Development Code contains a comprehensive lighting ordinance that appears similar to Flagstaff’s. The code addresses lighting controls, lumen caps, shielding and cut-off standards, sign illumination restrictions, and more. [17]

The code of ordinances for Clanton lacks a section devoted to lighting. Still, lights are mentioned in relation to other regulated items like signs and gas stations. The ordinance for gas stations includes a provision to prevent light trespass, but otherwise, outdoor lighting policy does not appear to be a priority for the city. [18]

Conclusions

Nighttime light emissions from the contiguous US show an increasing trend. Some states show drastic growth in radiance from 2014-2024, while others, particularly in the Midwest and Northeast, show more moderate changes.

After controlling for population and urban area, the nighttime lights of DarkSky-certified communities appear lower than those of similar cities that have not made outdoor lighting policy a priority. Bee Cave, the most recently designated Dark Sky place, is most similar to its paired community. It is possible that Bee Cave is still transitioning fixtures and practices to align with its newer commitment to dark sky preservation.

Limitations and further opportunities:

Radiance of nighttime lights is not a direct measure of light pollution. Some cities are switching from orange high pressure sodium lamps to white LEDs. VIIRS spectral range does not completely cover the blue light peak of many modern white LEDs, which can lead to somewhat misleading decreases in satellite-observed radiance [19].
Cities were paired based on 2023 ACS 5-year population estimates. As discussed, different rates of population growth should be considered. Further analysis could leverage more sophisticated use of population and land cover data. Census block groups could be rasterized to estimate the population within a nighttime lights grid cell.
One year of land cover data was assumed to represent each city’s land cover over the 10-year span, which is a simplification. Future analysis could incorporate a more nuanced use of land cover data (once AppEEARs processing is operational again).
Analysis could be expanded to additional cities.

References

[1] Chepesiuk, R. (2009). Missing the dark: health effects of light pollution. Environmental Health Perspectives, 117(1), A20+. https://www.jstor.org/stable/40066663.

[2] Hunter, T., (2023). The birth of DarkSky (IDA) and a lifelong mission fighting light pollution. DarkSky. https://darksky.org/news/the-birth-of-ida-and-a-lifelong-mission-fighting-light-pollution/.

[3] DarkSky International. (2025a). About DarkSky. https://darksky.org/about/#mission.

[4] DarkSky International (2025b). International Dark Sky Places. https://darksky.org/what-we-do/international-dark-sky-places/.

[5] Román, M. O., Wang, Z., Sun, Q., Kalb, V., Miller, S. D., Molthan, A., Schultz, L., Bell, J., Stokes, E. C., Pandey, B., Seto, K. C., Hall, D., Oda, T., Wolfe, R. E., Lin, G., Golpayegani, N., Devadiga, S., Davidson, C., Sarkar, S., Praderas, C., Schmaltz, J., Boller, R., Stevens, J., Ramos González, O. M., Padilla, E., Alonso, J., Detrés, Y., Armstrong, R., Miranda, I., Conte, Y., Marrero, N., MacManus, K., Esch, T., Masuoka, E. J. (2018). NASA’s Black Marble nighttime lights product suite. Remote Sensing of Environment, 210: 113-143, https://doi.org/10.1016/j.rse.2018.03.017.

[6] Marty, R., Vicente, G.S. (2025). Package ‘blackmarbler’. https://cran.r-project.org/web/packages/blackmarbler/blackmarbler.pdf.

[7] U.S. Census Bureau. 2024. 2024 Census TIGER/Line Shapefiles. https://data.census.gov

[8] Walker, K. , Rudis, B. (2025). Tigris: Load Census TIGER/Line shapefiles. R package version 2.2.1. https://cran.r-project.org/web/packages/tigris/index.html

[9] U.S. Census Bureau. (2023). Total Population. American Community Survey 5-Year Estimates Detailed Tables. Table B01003.

[10] Walker, K., Herman, M. (2025). tidycensus: Load US Census Boundary and Attribute Data as ‘tidyverse’ and ‘sf’-Ready Data Frames. R package version 1.7.3, https://walker-data.com/tidycensus/.

[11] Friedl, M., Sulla-Menashe, D. (2022). MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V061. NASA EOSDIS Land Processes Distributed Active Archive Center. Accessed 2025-11-29 from https://doi.org/10.5067/MODIS/MCD12Q1.061. Accessed November 29, 2025.

[12] Schultz, J. (2022). States shut out light pollution. National Conference of State Legislatures. https://www.ncsl.org/environment-and-natural-resources/states-shut-out-light-pollution

[13] Flagstaff, Arizona. (2011). Flagstaff Zoning Code, Division 10-50.70: Outdoor Lighting Standards. Accessed via DarkSky.org, https://darksky.org/places/flagstaff-arizona-dark-sky-community/.

[14] Scranton, Pennsylvania. (2023). City of Scranton Zoning Ordinance, Article 5, Section 5.10 - Lighting and Glare.

[15] Village of Homer Glen, Illinois. (2010). An Ordinance Regulating Outdoor Lighting in the Village of Homer Glen. Accessed via DarkSky.org, https://darksky.org/places/homer-glen-illinois-dark-sky-community/.

[16] Belvidere, Illinois. (2025). Belvidere, Illinois - Code of Ordinances, Chapter 150, Article 7, Section 150.707. - Exterior lighting standards. https://library.municode.com/il/belvidere/codes/code_of_ordinances?nodeId=BELVIDERE_MUNICIPAL_CODE_CH150ZO_ART7PEST_S150.707EXLIST

[17] City of Bee Cave, Texas. (2025). City of Bee Cave, TX Code, Title II, Chapter UDC, Article 6, Section 6.2 - Lighting. https://ecode360.com/40281623

[18] Clanton, Alabama. (2025). Code of Ordinances. https://library.municode.com/al/clanton/codes/code_of_ordinances?nodeId=12327

[19] Stare, Jurij; Kyba, Christopher (2025): Radiance Light Trends. GFZ Data Services. https://doi.org/10.5880/GFZ.1.4.2019.001.

Footnotes

Although land cover changes over the period from 2014 to 2024, for the purposes of this project, only one year of land cover data was considered for each community. Data availability is one reason for this simplification – annual data is not available for all study areas for all years, and the MODIS science team urges caution when using land cover layers for 2021 and beyond due to lack of updates to the classification training database [11]. Therefore, 2019 data was used for all cities except for Clanton, which used 2020 data.↩︎

--- title: "Starry Night?" subtitle: An exploration of nighttime lights and the Dark Sky Movement in the US author: Carly Staebell date: today date-format: long code-fold: true code-tools: toggle: true execute: message: false warning: false echo: true output: true --- # Introduction Electric lights have benefited society immensely since their invention in 1879. Modern work, recreation, and transportation are made possible largely due to artificial lighting. However, when artificial outdoor lighting becomes inefficient and unnecessary, it is known as light pollution, with negative impacts on both wildlife and humans [1]. Concerns about the impacts of light pollution on astronomical observation led to the 1988 formation of the nonprofit now known as DarkSky International [2]. Since then, DarkSky has grown into a recognized global authority on light pollution with a mission to safeguard communities and wildlife against light pollution [3]. The organization's International Dark Sky Places program allows communities, parks, and protected areas to seek certification for efforts to preserve and protect their dark skies [4]. While many communities implement outdoor lighting policies without seeking recognition, a DarkSky certification is an indicator of places that have made dark sky conservation a priority. This project aims to explore trends in nighttime light levels first from a national contiguous US perspective, and through more detailed case studies of individual cities. Three DarkSky certified communities [4] with different sizes and designation dates have been selected: * Flagstaff, Arizona: Population ~77,000, designated 2001; * Homer Glen, Illinois: Population ~25,000, designated 2011; * Bee Cave, Texas: Population ~9,000, designated 2023. These three communities were each paired with a community of similar size that has fewer, if any, ordinances controlling light pollution. Each community pair was examined to see if there are observable differences in nighttime radiance. # Materials and methods The following data were used. * NASA's Black Marble annual nighttime lights (VNP46A4 data product) [5] sourced via the `blackmarbler` package [6] + Spatial resolution: 500 m * US Census Bureau TIGER data [7] sourced via the `tigris` package [8] * US Census Bureau American Community Survey (ACS) population data [9] sourced via the `tidycensus` package [10] * MODIS Type 1 yearly land use/land cover (MCD12Q1 data product) sourced via AppEEARS [11] + Spatial resolution: 500 m ```{r load packages} #| output: false # Load libraries library(tidyverse) library(leaflet) library(kableExtra) library(htmlwidgets) library(widgetframe) knitr::opts_chunk$set(widgetframe_widgets_dir = 'widgets' ) knitr::opts_chunk$set(cache=TRUE) # cache the results for quick compiling library(blackmarbler) library(sf) library(terra) library(tigris) library(tidyterra) library(tidycensus) library(patchwork) library(piggyback) library(ggiraph) library(scales) ``` ## Download and clean all required data ### National and statewide nighttime lights A US Census TIGER shapefile was downloaded to provide geographic state boundaries. A separate script was used to download and save annual nighttime lights rasters for the contiguous US ('black_marble_data_download.R'). The script was also used to obtain extracted nighttime lights values aggregated at the national and state levels. Nighttime lights are provided as radiance values with units of $\frac{n W}{cm^2 sr}$. ```{r Download and process NTL data} #| output: false # Download and process CONUS NTL data # Define national spatial extent nonconus <- c("Guam", "Hawaii", "Alaska", "Commonwealth of the Northern Mariana Islands", "United States Virgin Islands", "American Samoa", "Puerto Rico") states_sf <- states() |> filter(!NAME %in% nonconus) |> select(NAME, geometry) |> st_transform(crs = "epsg:4326") # Dissolve individual state boundaries to get CONUS outline conus <- st_union(states_sf)|> st_as_sf() ## Read in national raster data for 2024 (to be used for leaflet map later) pb_download("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif", repo = "cstaebell/final-project-cstaebell") us_2024_rast <- rast("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif") |> mask(conus) # create log transformed version for plotting log_us_2024 <- app(us_2024_rast, fun = log1p) # Downsample log data for map (commented out since file written) #aggregate(log_us_2024, fact = 4, fun = "median", filename = "data/log_us_2024_downsampled.tif", overwrite = TRUE) #--- # Download and process annual US mean NTL values #--- years = 2014:2024 conus_means = numeric(length(years)) for (i in 1:length(years)) { pb_download(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t", years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell") annual_mean <- readRDS(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t", years[i], ".Rds", sep = "")) |> pull(ntl_mean) conus_means[i] <- annual_mean } conus_df <- data.frame(year = years, mean_ntl = conus_means) #--- # Download and process mean NTL values for states #--- state_means <- matrix(nrow = 49, ncol = length(years)) # Read in annual data for states and create a matrix of state means for (i in 1:length(years)) { pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t", years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell") annual_mean <- readRDS(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t", years[i], ".Rds", sep = "")) state_means[,i] <- annual_mean[,2] # check that state names are ordered the same way for each year if (i == 1) { state_names <- annual_mean[,1] } else { if (sum(state_names != annual_mean[,1]) != 0) { print("Different state names") } } } # Create data frame of state means and add percent change column state_means_df <- data.frame(state_means) colnames(state_means_df) <- paste("ntl_", years, sep = "") state_means_df <- state_means_df |> mutate(state = state_names, pct_change = (ntl_2024 - ntl_2014) / ntl_2014 * 100) # Combine state NTL data with geographic data states_ntl <- states_sf |> inner_join(state_means_df, by = c("NAME" = "state")) ``` ### Case Studies For case study analyses, population data and local geographies for each Dark Sky Place were downloaded from the US Census. ```{r Download census data} #| output: false # Download population data # US population to search for comparison cities us_pop <- get_acs(geography = "place", state = NULL, variables = "B01003_001", key = census_api_key) # Download local geographies: selected Dark Sky Places flagstaff_sf <- places(state = "AZ", cb = TRUE, year = 2024) |> filter(NAME == "Flagstaff") |> st_transform(crs = "epsg:4326") |> left_join(us_pop, by = "GEOID") homerglen_sf <- places(state = "IL", cb = TRUE, year = 2024) |> filter(NAME == "Homer Glen") |> st_transform(crs = "epsg:4326") |> left_join(us_pop, by = "GEOID") beecave_sf <- places(state = "TX", cb = TRUE, year = 2024) |> filter(NAME == "Bee Cave") |> st_transform(crs = "epsg:4326") |> left_join(us_pop, by = "GEOID") ``` ACS population data were filtered to return a listing of candidate municipalities with populations within 1% of each Dark Sky community population. Final Comparison cities were selected based on population match and state, with preference given to cities within states that are not known to have laws regarding outdoor lighting geared toward dark sky protection [12]. ```{r Process Census data} #| output: false # Get listing of suitable comparison cities flagstaff_sim_pop <- us_pop |> filter(estimate >= (flagstaff_sf$estimate - 0.005*flagstaff_sf$estimate) & estimate <= (flagstaff_sf$estimate + 0.005*flagstaff_sf$estimate)) |> separate(NAME, into = c("city", "state"), sep = ", ") homerglen_sim_pop <- us_pop |> filter(estimate >= (homerglen_sf$estimate - 0.005*homerglen_sf$estimate) & estimate <= (homerglen_sf$estimate + 0.005*homerglen_sf$estimate)) |> separate(NAME, into = c("city", "state"), sep = ", ") beecave_sim_pop <- us_pop |> filter(estimate >= (beecave_sf$estimate - 0.005*beecave_sf$estimate) & estimate <= (beecave_sf$estimate + 0.005*beecave_sf$estimate)) |> separate(NAME, into = c("city", "state"), sep = ", ") # Download/process comparison city data # Flagstaff counterpart scranton_sf <- places(state = "PA", cb = TRUE, year = 2024) |> filter(NAME == "Scranton") |> left_join(us_pop, by = "GEOID") |> st_transform(crs = "epsg:4326") # Homer Glen counterpart belvidere_sf <- places(state = "IL", cb = TRUE, year = 2024) |> filter(NAME == "Belvidere") |> left_join(us_pop, by = "GEOID") |> st_transform(crs = "epsg:4326") # Bee Cave counterpart clanton_sf <- places(state = "AL", cb = TRUE, year = 2024) |> filter(NAME == "Clanton") |> left_join(us_pop, by = "GEOID") |> st_transform(crs = "epsg:4326") city_names <- c("flagstaff", "homerglen", "beecave", "scranton", "belvidere", "clanton") ``` The final city pairs are as follows: | Pair | Community | Population (2023) | |------|-----------|------------| | A | Flagstaff, AZ | `r format(flagstaff_sf$estimate, big.mark = ",", scientific = FALSE)` | | A | Scranton, PA | `r format(scranton_sf$estimate, big.mark = ",", scientific = FALSE)` | | B | Homer Glen, IL | `r format(homerglen_sf$estimate, big.mark = ",", scientific = FALSE)` | | B | Belvidere, IL | `r format(belvidere_sf$estimate, big.mark = ",", scientific = FALSE)` | | C | Bee Cave, TX | `r format(beecave_sf$estimate, big.mark = ",", scientific = FALSE)` | | C | Clanton, AL | `r format(clanton_sf$estimate, big.mark = ",", scientific = FALSE)` | Annual population data for each city were downloaded from the ACS. The 5-year ACS was used for years 2014-2023. Since the 1-year ACS is only available for cities with populations > 65,000, 2024 data could only be obtained for Pair A cities. For the remaining cities, 2023 data was used for 2024. ```{r Annual population data} # 2014-2023 uses 5-year ACS population data annual_acs <- purrr::map(years[1:10], ~ get_acs(geography = "place", state = c("AZ", "IL", "TX", "PA", "AL"), variables = "B01003_001", key = census_api_key, year = .x, survey = "acs5")) |> map2(years[1:10], ~ mutate(.x, id = .y)) |> list_rbind() |> filter(GEOID %in% c(flagstaff_sf$GEOID, homerglen_sf$GEOID, beecave_sf$GEOID, scranton_sf$GEOID, belvidere_sf$GEOID, clanton_sf$GEOID)) |> separate(NAME, into = c("city", "state"), sep = ", ") |> rename(year = id, pop_est = estimate) |> mutate(city = case_when( city == "Flagstaff city" ~ "flagstaff", city == "Scranton city" ~ "scranton", city == "Homer Glen village" ~ "homerglen", city == "Belvidere city" ~ "belvidere", city == "Bee Cave city" ~ "beecave", city == "Clanton city" ~ "clanton", )) # 2024 requires 1-year ACS: only available for geographies with populations of 65000 and greater pop_2024 <- get_acs(geography = "place", state = c("AZ", "PA"), variables = "B01003_001", key = census_api_key, year = 2024, survey = "acs1") flagstaff_2024 <- pop_2024 |> filter(GEOID == flagstaff_sf$GEOID) |> separate(NAME, into = c("city", "state"), sep = ", ") |> rename(pop_est = estimate) |> mutate(year = 2024, city = "flagstaff") scranton_2024 <- pop_2024 |> filter(GEOID == scranton_sf$GEOID) |> separate(NAME, into = c("city", "state"), sep = ", ") |> rename(pop_est = estimate) |> mutate(year = 2024, city = "scranton") # Since data unavailable for smaller cities for 2024, impute 2023 data other_cities_2024 <- annual_acs |> filter(year == 2023 & city != "flagstaff" & city != "scranton") |> mutate(year = 2024) all_annual_acs <- annual_acs |> rbind(flagstaff_2024, scranton_2024, other_cities_2024) ``` After selecting cities, the annual US nighttime lights rasters were read in and subsetted to create a raster and dataframe of mean values for each city. Since simply looking at annual means for each city may be more of a proxy for overall development than for lighting policy, land use/land cover data[^1] were also brought in to help make comparisons across communities more meaningful. The land cover rasters were cropped to each city, and then the urban areas were isolated, allowing for calculation of annual mean nighttime lights from urban areas only. [^1]: Although land cover changes over the period from 2014 to 2024, for the purposes of this project, only one year of land cover data was considered for each community. Data availability is one reason for this simplification -- annual data is not available for all study areas for all years, and the MODIS science team urges caution when using land cover layers for 2021 and beyond due to lack of updates to the classification training database [11]. Therefore, 2019 data was used for all cities except for Clanton, which used 2020 data. ```{r Read in yearly national raster data} #| output: false # Nighttime lights and land cover data are in a single chunk to avoid external pointer errors when rendering #--- # Nighttime lights #--- # Subset nighttime lights rasters to each city city_sf_list <- list(flagstaff_sf, homerglen_sf, beecave_sf, scranton_sf, belvidere_sf, clanton_sf) city_ntl_rasters <- vector("list", length = length(city_names)) # Load files for each year for (i in 1:length(years)){ pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t", years[i], ".tif", sep = ""), repo = "cstaebell/final-project-cstaebell") rast_file <- paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t", years[i], ".tif", sep = "") year_rast <- rast(rast_file) # Loop over cities for (j in 1:length(city_sf_list)) { # Crop and mask to city extent current_rast <- year_rast |> terra::crop(city_sf_list[[j]]) |> mask(city_sf_list[[j]]) # Create list entry during first iteration, then append in subsequent iterations if (i == 1) { city_ntl_rasters[[j]] <- list(current_rast) } else { city_ntl_rasters[[j]][[length(city_ntl_rasters[[j]]) + 1]] <- current_rast } } } # Stack city lists city_ntl_rasters <- lapply(city_ntl_rasters, rast) # Create separate rasters for each city #city_ntl_var_names <- paste0(city_names, "_ntl_rast") #for (c in 1:length(city_ntl_var_names)) { # assign(city_ntl_var_names[c], city_ntl_rasters[[c]]) #} # Calculate annual means annual_ntl_means <- vector("list", length = length(city_names)) for (i in 1:length(city_ntl_rasters)) { ntl_vals <- data.frame(values(city_ntl_rasters[[i]])) |> summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |> pivot_longer(cols = t2014:t2024) |> mutate(name = years, city = city_names[i]) |> rename(year = name, mean_ntl = value) annual_ntl_means[[i]] <- ntl_vals } ntl_means_df <- list_rbind(annual_ntl_means) |> left_join(all_annual_acs, by = join_by(city == city, year == year)) |> mutate(percap_mean_ntl = mean_ntl / pop_est) #--- # Land cover data #--- # Download land cover data pb_download("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif", repo = "cstaebell/final-project-cstaebell") pb_download("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif", repo = "cstaebell/final-project-cstaebell") land_use_rast_2019 <- rast("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif") land_use_rast_2020 <- rast("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif") # Subset land cover rasters to each city land_use_rasters <- vector("list", length = length(city_sf_list)) urban_ntl_dfs <- vector("list", length = length(city_sf_list)) urban_ntl_vals <- vector("list", length = length(city_sf_list)) group_assign <- c("Dark Sky", "Dark Sky", "Dark Sky", "Comparison", "Comparison", "Comparison") pair_assign <- c("A", "B", "C", "A", "B", "C") for (i in 1:length(city_names)) { # Specify which data to use (2019 for all cities but Clanton) if (i == 6) { land_use_rast <- land_use_rast_2020 } else { land_use_rast <- land_use_rast_2019 } # Crop and mask to city extent current_rast <- land_use_rast |> terra::crop(city_sf_list[[i]]) |> mask(city_sf_list[[i]]) # Save to list for future use land_use_rasters[[i]] <- current_rast # Mask all but urban land use (land use code = 13) urban_area <- current_rast urban_area[urban_area != 13] <- NA city_urban <- city_ntl_rasters[[i]] |> mask(urban_area) # Calculate mean NTL for urban areas urban_ntl_dfs[[i]] <- data.frame(values(city_urban)) |> summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |> pivot_longer(cols = t2014:t2024) |> mutate(name = years, city = city_names[i]) |> rename(year = name, mean_urban_ntl = value) # Compile values for statistical tests urban_ntl_vals[[i]] <- data.frame(values(city_urban)) |> pivot_longer(cols = t2014:t2024) |> filter(!is.na(value)) |> mutate(city = city_names[[i]], group = group_assign[i], pair = pair_assign[i]) } # Combine all mean NTL data into single data frame all_ntl_means_df <- ntl_means_df |> left_join(list_rbind(urban_ntl_dfs), join_by(year == year, city == city)) |> mutate(mean_percap_urban_ntl = mean_urban_ntl / pop_est, pair = case_when( city == "flagstaff" | city == "scranton" ~ "A", city == "homerglen" | city == "belvidere" ~ "B", city == "beecave" | city == "clanton" ~ "C"), group = case_when( city == "flagstaff" | city == "homerglen" | city == "beecave" ~ "Dark Sky", city == "scranton" | city == "belvidere" | city == "clanton" ~ "Comparison") ) # Prepare data for sample land cover plots # Define land cover encoding Land_Cover_Type_1 <- c( "Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2, "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5, "Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9, "Grasslands" = 10, "Permanent wetlands" = 11, "Croplands" = 12, "Urban & built-up" = 13, "Cropland/Natural vegetation mosaic" = 14, "Snow & ice" = 15, "Barren/Sparsely vegetated" = 16 ) lcd <- data.frame( ID = Land_Cover_Type_1, landcover = names(Land_Cover_Type_1), stringsAsFactors = FALSE) # Create data frames for plotting flagstaff_lulc <- as.data.frame(land_use_rasters[[1]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |> mutate(city = "Flagstaff") scranton_lulc <- as.data.frame(land_use_rasters[[4]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |> mutate(city = "Scranton") homerglen_lulc <- as.data.frame(land_use_rasters[[2]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |> mutate(city = "Homer Glen") belvidere_lulc <- as.data.frame(land_use_rasters[[5]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |> mutate(city = "Belvidere") beecave_lulc <- as.data.frame(land_use_rasters[[3]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |> mutate(city = "Bee Cave") clanton_lulc <- as.data.frame(land_use_rasters[[6]], xy = TRUE) |> left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2020001000000_aid0001 == ID)) |> rename(lc_code = MCD12Q1.061_LC_Type1_doy2020001000000_aid0001) |> mutate(city = "Clanton") all_lulc <- bind_rows(flagstaff_lulc, scranton_lulc, homerglen_lulc, belvidere_lulc, beecave_lulc, clanton_lulc) all_lulc$city <- factor(all_lulc$city, levels = c("Flagstaff", "Scranton", "Homer Glen", "Belvidere", "Bee Cave", "Clanton")) ``` Statistical testing was done on the nonaggregated urban lights values. Although there are thousands of samples, radiance data exhibit a heavy right skew. The Shapiro-Wilk test was used to confirm non-normality, and since the data are not normal, a Wilcoxon test was used for each community pair to test if the Comparison group radiance is greater than that of the Dark Sky group. ```{r Stat test} #| results: hide # Prep values for statistical testing pair_a_uvals <- bind_rows(urban_ntl_vals[[1]], urban_ntl_vals[[4]]) pair_b_uvals <- bind_rows(urban_ntl_vals[[2]], urban_ntl_vals[[5]]) pair_c_uvals <- bind_rows(urban_ntl_vals[[3]], urban_ntl_vals[[6]]) for (i in 1:length(urban_ntl_vals)) { shapiro_test <- shapiro.test(urban_ntl_vals[[i]]$value) if (shapiro_test$p.value < 0.05) { print("Reject normality") } else { print("Fail to reject normality") } } wilcox_a <- wilcox.test(value ~ group, data = pair_a_uvals, alternative = "greater") wilcox_b <- wilcox.test(value ~ group, data = pair_b_uvals, alternative = "greater") wilcox_c <- wilcox.test(value ~ group, data = pair_c_uvals, alternative = "greater") # Repeat analysis for only 2021 - 2024 due to population trends pair_a_recent <- pair_a_uvals |> filter(name %in% c("t2021", "t2022", "t2023", "t2024")) pair_b_recent <- pair_b_uvals |> filter(name %in% c("t2021", "t2022", "t2023", "t2024")) pair_c_recent <- pair_c_uvals |> filter(name %in% c("t2021", "t2022", "t2023", "t2024")) wilcox_a_recent <- wilcox.test(value ~ group, data = pair_a_recent, alternative = "greater") wilcox_b_recent <- wilcox.test(value ~ group, data = pair_b_recent, alternative = "greater") wilcox_c_recent <- wilcox.test(value ~ group, data = pair_c_recent, alternative = "greater") ``` # Results ## National Trends Mean nighttime light radiance for the entire contiguous US is plotted below, where a generally increasing trend can be observed over the past ten years. ```{r Plot CONUS mean radiance vs year} ggplot(conus_df) + geom_point(aes(x = year, y = mean_ntl)) + geom_smooth(aes(x = year, y = mean_ntl)) + labs(x = "Year", y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"), title = "Annual US Mean Radiance (2014 - 2024)") + theme_light() ``` State means for 2024 are shown in the figure below. Annual radiance tends to be higher for eastern states, which reflects known patterns of population density. ```{r national mean choropleth} # Create ggplot map for ggiraph st_2024_mean_map <- ggplot() + geom_sf_interactive(data = filter(states_ntl, NAME != "District of Columbia"), aes(fill = ntl_2024, tooltip = paste(NAME, "\n", round(ntl_2024, 2)), data_id = NAME)) + theme_void() + theme(text = element_text(family = "Arial")) + scale_fill_viridis_c() + labs(title = "Mean Nighttime Lights by State, 2024", fill = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"), caption = "Data: NASA Black Marble & US Census Bureau") # Transform into interactive ggiraph object girafe(ggobj = st_2024_mean_map, options = list(opts_hover(css = "stroke:cyan;stroke-width:2px"), opts_hover_inv(css = "opacity:0.8;"), opts_tooltip(css = "background-color:white;color:black;padding:15px; border-radius:10px; font-family:Arial;font-size:15px;" ))) ``` In the period from 2014 to 2024, most states experienced an increase in detected nighttime radiance. Texas, Montana, Idaho, and Utah all show particularly large increases upwards of 40%. More research is required to explain these trends, but possible drivers include recent population growth and urban development. On the other hand, several states have stayed fairly constant or even decreased over the 10-year span. ```{r state pct change} st_pct_change_map <- ggplot() + geom_sf_interactive(data = filter(states_ntl, NAME != "District of Columbia"), aes(fill = pct_change, tooltip = paste(NAME, "\n", round(pct_change, 1),"%", "\n2014: ", round(ntl_2014, 2), "\n2024: ", round(ntl_2024, 2)), data_id = NAME)) + theme_void() + theme(text = element_text(family = "Arial")) + scale_fill_gradient2(low = "navy", mid = "white", high = "darkred", midpoint = 0) + labs(title = "Mean Radiance Percent Change (2014 - 2024)", fill = "% Change", caption = "Data: NASA Black Marble & US Census Bureau") girafe(ggobj = st_pct_change_map, options = list(opts_hover(css = "stroke:black;stroke-width:2px"), opts_hover_inv(css = "opacity:0.8;"), opts_tooltip(css = "background-color:white;color:black;padding:15px; border-radius:10px; font-family:Arial;font-size:15px;"))) ``` ## Case Studies ```{r leaflet map} # Leaflet prep fs_coord <- st_centroid(flagstaff_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Flagstaff", pair = "A", group = "Dark Sky") sc_coord <- st_centroid(scranton_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Scranton", pair = "A", group = "Comparison") hg_coord <- st_centroid(homerglen_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Homer Glen", pair = "B", group = "Dark Sky") bv_coord <- st_centroid(belvidere_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Belvidere", pair = "B", group = "Comparison") bc_coord <- st_centroid(beecave_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Bee Cave", pair = "C", group = "Dark Sky") cl_coord <- st_centroid(clanton_sf) |> st_coordinates() |> data.frame() |> mutate(city = "Clanton", pair = "C", group = "Comparison") markers_df <- bind_rows(fs_coord, sc_coord, hg_coord, bv_coord, bc_coord, cl_coord) markers_df$pair <- as.factor(markers_df$pair) log_us_2024_ds <- rast("data/log_us_2024_downsampled.tif") vals <- values(log_us_2024_ds, mat = FALSE) midpoint <- 3.5 breaks <- c(min(vals, na.rm = TRUE), seq(min(vals, na.rm = TRUE), midpoint, length.out = 5)[-1], seq(midpoint, max(vals, na.rm = TRUE), length.out = 5)[-1]) pal <- colorBin(palette = c("black", "yellow", "red"), domain = vals, bins = breaks, na.color = "transparent") city_color_pal <- colorFactor(palette = c("firebrick", "royalblue3", "green4"), domain = markers_df$pair) city_popup <- paste0("", markers_df$city, " ", "", markers_df$group, " Community ", "Pair: ", markers_df$pair) # Leaflet map leaflet() |> addProviderTiles("CartoDB") |> addRasterImage(log_us_2024_ds, group = "US NTL", colors = pal) |> addCircleMarkers(data = markers_df, ~X, ~Y, fillColor = ~city_color_pal(pair), fillOpacity = 0.9, color = "black", label = ~as.character(city), popup = city_popup) ``` From the plot below, four of the six case study cities exhibit an increasing mean radiance trend similar to that of the national mean. However, the Pair A cities, Flagstaff and Scranton, show decreasing tendencies. The large difference in radiance values within pairs is another notable feature of this plot. ```{r compare city annual means} # Prep data frame so that plots show in desired order all_ntl_means_df$pair <- as.factor(all_ntl_means_df$pair) pair_names <- c("A" = "Flagstaff & Scranton", "B" = "Homer Glen & Belvidere", "C" = "Bee Cave & Clanton") # Plot of annual means ggplot(data = all_ntl_means_df) + geom_line(aes(x = year, y = mean_ntl, color = pair, linetype = group), linewidth = 1) + facet_wrap(~ pair, labeller = labeller(pair = pair_names)) + scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) + scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) + labs(x = "Year", y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"), title = "Mean Annual Radiance", linetype = "City Type", color = "Pair") + theme_light() + theme(strip.text.x = element_text(color = "black"), plot.title = element_text(hjust = 0.5), aspect.ratio = 1.8) ``` Simply looking at the annual means provides interesting information, but fails to take into account two factors that may influence a city's observed nighttime light emissions: population over time and land cover. Annual population data reveals that Flagstaff and Bee Cave have both experienced steep population growth from 2014 to 2024. While Belvidere shows a generally decreasing trend, both Belvidere and Homer Glen had relatively stable populations with less than 5% change. Scranton and Clanton's populations have also remained quite stable. ```{r Population trends} ggplot(data = all_ntl_means_df) + geom_line(aes(x = year, y = pop_est, color = pair, linetype = group), linewidth = 1) + facet_wrap(~ pair, labeller = labeller(pair = pair_names), scale = "free_y") + scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) + scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) + labs(x = "Year", y = "Population", title = "Annual Populations", linetype = "City Type", color = "Pair") + theme_light() + theme(strip.text.x = element_text(color = "black"), plot.title = element_text(hjust = 0.5), axis.text.y = element_text(size = 8), aspect.ratio = 1.8) ``` Land cover differences between paired cities helps to explain some of the features and differences observed in the plot of mean nighttime lights versus time. Flagstaff contains much more undeveloped land like grassland and savannas, whereas most of Scranton is classified as urban and built up, which likely accounts for the much lower mean values for Flagstaff. Similarly, Clanton has a very small amount of land classified as urban & built-up, which probably contributes to a low mean radiance. ```{r Land cover plots} ggplot(all_lulc, aes(x = x, y = y)) + geom_tile(aes(fill = landcover)) + facet_wrap(~ city, nrow = 3, scales = "free") + coord_quickmap() + theme_light() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), strip.text.x = element_text(color = "black")) + labs(fill = "Land Cover", caption = "Data: NASA & US Census Bureau") + scale_fill_viridis_d() ``` Given the imbalances in land cover between different cities, it is useful for the sake of comparison to restrict analysis to urban areas. Proceeding this way introduces the assumption that any lighting policy effects outside of non-built-up areas are dwarfed by those on urban areas. From the plots of urban mean radiance, we can see that the radiance for the Dark Sky Places is generally lower than for their paired cities. ```{r Urban NTL comparison} # Urban mean NTL plot ggplot(data = all_ntl_means_df) + geom_line(aes(x = year, y = mean_urban_ntl, color = pair, linetype = group), linewidth = 1) + facet_wrap(~ pair, labeller = labeller(pair = pair_names), scales = "free_y") + scale_linetype_manual(values = c("Dark Sky" = "solid", "Comparison" = "dashed")) + scale_color_manual(values = c("A" = "firebrick", "B" = "royalblue3", "C" = "green4")) + labs(x = "Year", y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"), title = "Mean Annual Radiance from Urban Areas", linetype = "City Type", color = "Pair") + theme_light() + theme(strip.text.x = element_text(color = "black"), plot.title = element_text(hjust = 0.5), aspect.ratio = 1.8) ``` The results of the Wilcoxon test on the urban radiance values suggest that the observed differences between the Dark Sky communities and their Comparison counterparts are statistically significant. Given the population discrepancies between the Pair A and Pair C cities in the earlier analysis years, analysis was done on the entire set of data and on a subset of data from 2021-2024, when the populations became more comparable. | Pair | p-value (2014-2024) | p-value (2021-2024) | |------|---------------------|---------------------| | A | `r pvalue(wilcox_a$p.value, accuracy = 0.0001)` | `r pvalue(wilcox_a_recent$p.value, accuracy = 0.0001)` | | B | `r pvalue(wilcox_b$p.value, accuracy = 0.0001)` | `r pvalue(wilcox_b_recent$p.value, accuracy = 0.0001)` | | C | `r round(wilcox_c$p.value, digits = 5)` | `r round(wilcox_c_recent$p.value, digits = 5)` | ### Discussion of lighting policies #### Pair A As the home of two major observatories, Flagstaff has been at the forefront of outdoor lighting policy since 1958. The city's lighting ordinance aims to encourage lighting practices and systems that will minimize light pollution and trespass while still maintaining nighttime safety, security, and productivity. The code includes comprehensive standards for maximum total outdoor light output, shielding requirements, time limits, sign illumination, and lamp types. [13] Scranton's Zoning Ordinance includes one short section devoted to lighting, which mainly focuses on curbing light trespass and nuisance glare. Maximum illumination levels are defined for two lighting levels. [14] #### Pair B Homer Glen's outdoor lighting ordinance is rooted in the desire to preserve rural character and aesthetic value. The ordinance establishes maximum light output and inclination for various outdoor lighting applications, permitted hours for lighting, shielding requirements, and bulb types. [15] Similar to Scranton, Belvidere includes exterior lighting standards in their code of ordinances but does so for the purposes of promoting traffic safety and preventing nuisances. The ordinance sets limits on illumination intensity and hours of illumination outside of buildings, plus some requirements for light orientation and shielding. However, these provisions are not as strict or clearly defined as those in the Homer Glen ordinanace. [16] #### Pair C Bee Cave's Unified Development Code contains a comprehensive lighting ordinance that appears similar to Flagstaff's. The code addresses lighting controls, lumen caps, shielding and cut-off standards, sign illumination restrictions, and more. [17] The code of ordinances for Clanton lacks a section devoted to lighting. Still, lights are mentioned in relation to other regulated items like signs and gas stations. The ordinance for gas stations includes a provision to prevent light trespass, but otherwise, outdoor lighting policy does not appear to be a priority for the city. [18] # Conclusions Nighttime light emissions from the contiguous US show an increasing trend. Some states show drastic growth in radiance from 2014-2024, while others, particularly in the Midwest and Northeast, show more moderate changes. After controlling for population and urban area, the nighttime lights of DarkSky-certified communities appear lower than those of similar cities that have not made outdoor lighting policy a priority. Bee Cave, the most recently designated Dark Sky place, is most similar to its paired community. It is possible that Bee Cave is still transitioning fixtures and practices to align with its newer commitment to dark sky preservation. ### Limitations and further opportunities: * Radiance of nighttime lights is not a direct measure of light pollution. Some cities are switching from orange high pressure sodium lamps to white LEDs. VIIRS spectral range does not completely cover the blue light peak of many modern white LEDs, which can lead to somewhat misleading decreases in satellite-observed radiance [19]. * Cities were paired based on 2023 ACS 5-year population estimates. As discussed, different rates of population growth should be considered. Further analysis could leverage more sophisticated use of population and land cover data. Census block groups could be rasterized to estimate the population within a nighttime lights grid cell. * One year of land cover data was assumed to represent each city's land cover over the 10-year span, which is a simplification. Future analysis could incorporate a more nuanced use of land cover data (once AppEEARs processing is operational again). * Analysis could be expanded to additional cities. # References [1] Chepesiuk, R. (2009). Missing the dark: health effects of light pollution. Environmental Health Perspectives, 117(1), A20+. <https://www.jstor.org/stable/40066663>. [2] Hunter, T., (2023). The birth of DarkSky (IDA) and a lifelong mission fighting light pollution. DarkSky. <https://darksky.org/news/the-birth-of-ida-and-a-lifelong-mission-fighting-light-pollution/>. [3] DarkSky International. (2025a). About DarkSky. <https://darksky.org/about/#mission>. [4] DarkSky International (2025b). International Dark Sky Places. <https://darksky.org/what-we-do/international-dark-sky-places/>. [5] Román, M. O., Wang, Z., Sun, Q., Kalb, V., Miller, S. D., Molthan, A., Schultz, L., Bell, J., Stokes, E. C., Pandey, B., Seto, K. C., Hall, D., Oda, T., Wolfe, R. E., Lin, G., Golpayegani, N., Devadiga, S., Davidson, C., Sarkar, S., Praderas, C., Schmaltz, J., Boller, R., Stevens, J., Ramos González, O. M., Padilla, E., Alonso, J., Detrés, Y., Armstrong, R., Miranda, I., Conte, Y., Marrero, N., MacManus, K., Esch, T., Masuoka, E. J. (2018). NASA's Black Marble nighttime lights product suite. Remote Sensing of Environment, 210: 113-143, <https://doi.org/10.1016/j.rse.2018.03.017>. [6] Marty, R., Vicente, G.S. (2025). Package 'blackmarbler'. <https://cran.r-project.org/web/packages/blackmarbler/blackmarbler.pdf>. [7] U.S. Census Bureau. 2024. 2024 Census TIGER/Line Shapefiles. <https://data.census.gov> [8] Walker, K. , Rudis, B. (2025). Tigris: Load Census TIGER/Line shapefiles. R package version 2.2.1. <https://cran.r-project.org/web/packages/tigris/index.html> [9] U.S. Census Bureau. (2023). Total Population. American Community Survey 5-Year Estimates Detailed Tables. Table B01003. [10] Walker, K., Herman, M. (2025). tidycensus: Load US Census Boundary and Attribute Data as 'tidyverse' and 'sf'-Ready Data Frames. R package version 1.7.3, <https://walker-data.com/tidycensus/>. [11] Friedl, M., Sulla-Menashe, D. (2022). MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V061. NASA EOSDIS Land Processes Distributed Active Archive Center. Accessed 2025-11-29 from <https://doi.org/10.5067/MODIS/MCD12Q1.061>. Accessed November 29, 2025. [12] Schultz, J. (2022). States shut out light pollution. National Conference of State Legislatures. <https://www.ncsl.org/environment-and-natural-resources/states-shut-out-light-pollution> [13] Flagstaff, Arizona. (2011). Flagstaff Zoning Code, Division 10-50.70: Outdoor Lighting Standards. Accessed via DarkSky.org, <https://darksky.org/places/flagstaff-arizona-dark-sky-community/>. [14] Scranton, Pennsylvania. (2023). City of Scranton Zoning Ordinance, Article 5, Section 5.10 - Lighting and Glare. [15] Village of Homer Glen, Illinois. (2010). An Ordinance Regulating Outdoor Lighting in the Village of Homer Glen. Accessed via DarkSky.org, <https://darksky.org/places/homer-glen-illinois-dark-sky-community/>. [16] Belvidere, Illinois. (2025). Belvidere, Illinois - Code of Ordinances, Chapter 150, Article 7, Section 150.707. - Exterior lighting standards. <https://library.municode.com/il/belvidere/codes/code_of_ordinances?nodeId=BELVIDERE_MUNICIPAL_CODE_CH150ZO_ART7PEST_S150.707EXLIST> [17] City of Bee Cave, Texas. (2025). City of Bee Cave, TX Code, Title II, Chapter UDC, Article 6, Section 6.2 - Lighting. <https://ecode360.com/40281623> [18] Clanton, Alabama. (2025). Code of Ordinances. <https://library.municode.com/al/clanton/codes/code_of_ordinances?nodeId=12327> [19] Stare, Jurij; Kyba, Christopher (2025): Radiance Light Trends. GFZ Data Services. <https://doi.org/10.5880/GFZ.1.4.2019.001>.